0510Lesson

library(sf)

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(raster)

## Loading required package: sp

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(stringr) # for working with strings (pattern matching)
library(tidyr)   # for unite() and separate()

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:raster':
## 
##     extract

library(spData)

world

## Simple feature collection with 177 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 11
##    iso_a2 name_long continent region_un subregion type  area_km2     pop lifeExp
##    <chr>  <chr>     <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
##  1 FJ     Fiji      Oceania   Oceania   Melanesia Sove…   1.93e4  8.86e5    70.0
##  2 TZ     Tanzania  Africa    Africa    Eastern … Sove…   9.33e5  5.22e7    64.2
##  3 EH     Western … Africa    Africa    Northern… Inde…   9.63e4 NA         NA  
##  4 CA     Canada    North Am… Americas  Northern… Sove…   1.00e7  3.55e7    82.0
##  5 US     United S… North Am… Americas  Northern… Coun…   9.51e6  3.19e8    78.8
##  6 KZ     Kazakhst… Asia      Asia      Central … Sove…   2.73e6  1.73e7    71.6
##  7 UZ     Uzbekist… Asia      Asia      Central … Sove…   4.61e5  3.08e7    71.0
##  8 PG     Papua Ne… Oceania   Oceania   Melanesia Sove…   4.65e5  7.76e6    65.2
##  9 ID     Indonesia Asia      Asia      South-Ea… Sove…   1.82e6  2.55e8    68.9
## 10 AR     Argentina South Am… Americas  South Am… Sove…   2.78e6  4.30e7    76.3
## # … with 167 more rows, and 2 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>

coffee_data

## # A tibble: 47 x 3
##    name_long                coffee_production_2016 coffee_production_2017
##    <chr>                                     <int>                  <int>
##  1 Angola                                       NA                     NA
##  2 Bolivia                                       3                      4
##  3 Brazil                                     3277                   2786
##  4 Burundi                                      37                     38
##  5 Cameroon                                      8                      6
##  6 Central African Republic                     NA                     NA
##  7 Congo, Dem. Rep. of                           4                     12
##  8 Colombia                                   1330                   1169
##  9 Costa Rica                                   28                     32
## 10 Côte d'Ivoire                               114                    130
## # … with 37 more rows

world_coffee = left_join(world, coffee_data, by = "name_long")  #2つのdata.frameの共通点を元にくっつける。byで明示
#> Joining, by = "name_long"
class(world_coffee)

## [1] "sf"         "tbl_df"     "tbl"        "data.frame"

#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"
world_coffee

## Simple feature collection with 177 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 13
##    iso_a2 name_long continent region_un subregion type  area_km2     pop lifeExp
##    <chr>  <chr>     <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
##  1 FJ     Fiji      Oceania   Oceania   Melanesia Sove…   1.93e4  8.86e5    70.0
##  2 TZ     Tanzania  Africa    Africa    Eastern … Sove…   9.33e5  5.22e7    64.2
##  3 EH     Western … Africa    Africa    Northern… Inde…   9.63e4 NA         NA  
##  4 CA     Canada    North Am… Americas  Northern… Sove…   1.00e7  3.55e7    82.0
##  5 US     United S… North Am… Americas  Northern… Coun…   9.51e6  3.19e8    78.8
##  6 KZ     Kazakhst… Asia      Asia      Central … Sove…   2.73e6  1.73e7    71.6
##  7 UZ     Uzbekist… Asia      Asia      Central … Sove…   4.61e5  3.08e7    71.0
##  8 PG     Papua Ne… Oceania   Oceania   Melanesia Sove…   4.65e5  7.76e6    65.2
##  9 ID     Indonesia Asia      Asia      South-Ea… Sove…   1.82e6  2.55e8    68.9
## 10 AR     Argentina South Am… Americas  South Am… Sove…   2.78e6  4.30e7    76.3
## # … with 167 more rows, and 4 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>, coffee_production_2016 <int>,
## #   coffee_production_2017 <int>

plot(world_coffee["coffee_production_2017"])

world_coffee_inner <- inner_join(world, coffee_data)  #inner_join -> left_joinでNAになった要素を弾く

## Joining, by = "name_long"

world_coffee_inner

## Simple feature collection with 45 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -117.1278 ymin: -33.76838 xmax: 156.02 ymax: 35.49401
## Geodetic CRS:  WGS 84
## # A tibble: 45 x 13
##    iso_a2 name_long  continent region_un subregion type  area_km2    pop lifeExp
##    <chr>  <chr>      <chr>     <chr>     <chr>     <chr>    <dbl>  <dbl>   <dbl>
##  1 TZ     Tanzania   Africa    Africa    Eastern … Sove…  932746. 5.22e7    64.2
##  2 PG     Papua New… Oceania   Oceania   Melanesia Sove…  464520. 7.76e6    65.2
##  3 ID     Indonesia  Asia      Asia      South-Ea… Sove… 1819251. 2.55e8    68.9
##  4 KE     Kenya      Africa    Africa    Eastern … Sove…  590837. 4.60e7    66.2
##  5 DO     Dominican… North Am… Americas  Caribbean Sove…   48158. 1.04e7    73.5
##  6 TL     Timor-Les… Asia      Asia      South-Ea… Sove…   14715. 1.21e6    68.3
##  7 MX     Mexico     North Am… Americas  Central … Sove… 1969480. 1.24e8    76.8
##  8 BR     Brazil     South Am… Americas  South Am… Sove… 8508557. 2.04e8    75.0
##  9 BO     Bolivia    South Am… Americas  South Am… Sove… 1085270. 1.06e7    68.4
## 10 PE     Peru       South Am… Americas  South Am… Sove… 1309700. 3.10e7    74.5
## # … with 35 more rows, and 4 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>, coffee_production_2016 <int>,
## #   coffee_production_2017 <int>

plot(world_coffee_inner["coffee_production_2017"])

world_new = world # do not overwrite our original data
world_new$pop_dens = world_new$pop / world_new$area_km2 #計算されたpop_densをworld_newに追加
world_new

## Simple feature collection with 177 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 12
##    iso_a2 name_long continent region_un subregion type  area_km2     pop lifeExp
##  * <chr>  <chr>     <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
##  1 FJ     Fiji      Oceania   Oceania   Melanesia Sove…   1.93e4  8.86e5    70.0
##  2 TZ     Tanzania  Africa    Africa    Eastern … Sove…   9.33e5  5.22e7    64.2
##  3 EH     Western … Africa    Africa    Northern… Inde…   9.63e4 NA         NA  
##  4 CA     Canada    North Am… Americas  Northern… Sove…   1.00e7  3.55e7    82.0
##  5 US     United S… North Am… Americas  Northern… Coun…   9.51e6  3.19e8    78.8
##  6 KZ     Kazakhst… Asia      Asia      Central … Sove…   2.73e6  1.73e7    71.6
##  7 UZ     Uzbekist… Asia      Asia      Central … Sove…   4.61e5  3.08e7    71.0
##  8 PG     Papua Ne… Oceania   Oceania   Melanesia Sove…   4.65e5  7.76e6    65.2
##  9 ID     Indonesia Asia      Asia      South-Ea… Sove…   1.82e6  2.55e8    68.9
## 10 AR     Argentina South Am… Americas  South Am… Sove…   2.78e6  4.30e7    76.3
## # … with 167 more rows, and 3 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>, pop_dens <dbl>

world_new2 <- world %>% 
  mutate(pop_dens = pop/area_km2)
world_new2  #上と同じ

## Simple feature collection with 177 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 12
##    iso_a2 name_long continent region_un subregion type  area_km2     pop lifeExp
##  * <chr>  <chr>     <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
##  1 FJ     Fiji      Oceania   Oceania   Melanesia Sove…   1.93e4  8.86e5    70.0
##  2 TZ     Tanzania  Africa    Africa    Eastern … Sove…   9.33e5  5.22e7    64.2
##  3 EH     Western … Africa    Africa    Northern… Inde…   9.63e4 NA         NA  
##  4 CA     Canada    North Am… Americas  Northern… Sove…   1.00e7  3.55e7    82.0
##  5 US     United S… North Am… Americas  Northern… Coun…   9.51e6  3.19e8    78.8
##  6 KZ     Kazakhst… Asia      Asia      Central … Sove…   2.73e6  1.73e7    71.6
##  7 UZ     Uzbekist… Asia      Asia      Central … Sove…   4.61e5  3.08e7    71.0
##  8 PG     Papua Ne… Oceania   Oceania   Melanesia Sove…   4.65e5  7.76e6    65.2
##  9 ID     Indonesia Asia      Asia      South-Ea… Sove…   1.82e6  2.55e8    68.9
## 10 AR     Argentina South Am… Americas  South Am… Sove…   2.78e6  4.30e7    76.3
## # … with 167 more rows, and 3 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>, pop_dens <dbl>

world %>% transmute(pop_dens = pop/area_km2)  #計算されたものだけ残る

## Simple feature collection with 177 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 2
##    pop_dens                                                                 geom
##  *    <dbl>                                                   <MULTIPOLYGON [°]>
##  1    45.9  (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135, 178.7251 -17.0…
##  2    56.0  (((33.90371 -0.95, 34.07262 -1.05982, 37.69869 -3.09699, 37.7669 -3…
##  3    NA    (((-8.66559 27.65643, -8.665124 27.58948, -8.6844 27.39574, -8.6872…
##  4     3.54 (((-122.84 49, -122.9742 49.00254, -124.9102 49.98456, -125.6246 50…
##  5    33.5  (((-122.84 49, -120 49, -117.0312 49, -116.0482 49, -113 49, -110.0…
##  6     6.33 (((87.35997 49.21498, 86.59878 48.54918, 85.76823 48.45575, 85.7204…
##  7    66.7  (((55.96819 41.30864, 55.92892 44.99586, 58.50313 45.5868, 58.68999…
##  8    16.7  (((141.0002 -2.600151, 142.7352 -3.289153, 144.584 -3.861418, 145.2…
##  9   140.   (((141.0002 -2.600151, 141.0171 -5.859022, 141.0339 -9.117893, 140.…
## 10    15.4  (((-68.63401 -52.63637, -68.25 -53.1, -67.75 -53.85, -66.45 -54.45,…
## # … with 167 more rows

world %>% transmute(pop_dens = pop/area_km2) %>% st_drop_geometry()

## # A tibble: 177 x 1
##    pop_dens
##  *    <dbl>
##  1    45.9 
##  2    56.0 
##  3    NA   
##  4     3.54
##  5    33.5 
##  6     6.33
##  7    66.7 
##  8    16.7 
##  9   140.  
## 10    15.4 
## # … with 167 more rows

#raster画像を自分で作成
elev = raster(nrows = 6, ncols = 6, res = 0.5,
              xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
              vals = 1:36)
plot(elev)

grain_order = c("clay", "silt", "sand")
grain_char = sample(grain_order, 36, replace = TRUE)
grain_fact = factor(grain_char, levels = grain_order)
grain = raster(nrows = 6, ncols = 6, res = 0.5, 
               xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
               vals = grain_fact)
plot(grain)

# row 1, column 1
elev[1, 1]

## [1] 1

# cell ID 1
elev[1]

## [1] 1

r_stack = stack(elev, grain)
r_stack

## class      : RasterStack 
## dimensions : 6, 6, 36, 2  (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : layer.1, layer.2 
## min values :       1,       1 
## max values :      36,       3

names(r_stack) = c("elev", "grain")
# three ways to extract a layer of a stack
raster::subset(r_stack, "elev")

## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : elev 
## values     : 1, 36  (min, max)

r_stack[["elev"]]

## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : elev 
## values     : 1, 36  (min, max)

r_stack$elev

## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : elev 
## values     : 1, 36  (min, max)

elev[1, 1] = 0  #値の代入
elev[]

##  [1]  0  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36

cellStats(elev, sd)

## [1] 10.58432

summary(elev)

##         layer
## Min.     0.00
## 1st Qu.  9.75
## Median  18.50
## 3rd Qu. 27.25
## Max.    36.00
## NA's     0.00

library(sf)
library(raster)
library(dplyr)
library(spData)

canterbury = nz %>% filter(Name == "Canterbury")
canterbury_height = nz_height[canterbury, ]

sel_sgbp = st_intersects(x = nz_height, y = canterbury,sparse = FALSE) #sparseでTRUE or FALSEで表示
class(sel_sgbp)

## [1] "matrix" "array"

#> [1] "sgbp" "list"
sel_logical = lengths(sel_sgbp) > 0
canterbury_height2 = nz_height[sel_logical, ]

# create a polygon
a_poly = st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1))))
a = st_sfc(a_poly)
# create a line
l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 2))
l = st_sfc(l_line)
# create points
p_matrix = matrix(c(0.5, 1, -1, 0, 0, 1, 0.5, 1), ncol = 2)
p_multi = st_multipoint(x = p_matrix)
p = st_cast(st_sfc(p_multi), "POINT")

sel = st_is_within_distance(p, a, dist = 0.9) # can only return a sparse matrix
lengths(sel) > 0

## [1]  TRUE  TRUE FALSE  TRUE

#> [1]  TRUE  TRUE FALSE  TRUE

set.seed(2018) # set seed for reproducibility
(bb_world = st_bbox(world)) # the world's bounds

##       xmin       ymin       xmax       ymax 
## -180.00000  -90.00000  180.00000   83.64513

#>   xmin   ymin   xmax   ymax 
#> -180.0  -90.0  180.0   83.6
random_df = tibble(
  x = runif(n = 10, min = bb_world[1], max = bb_world[3]),
  y = runif(n = 10, min = bb_world[2], max = bb_world[4])
)
random_points = random_df %>% 
  st_as_sf(coords = c("x", "y")) %>% # set coordinates
  st_set_crs(4326) # set geographic CRS

world_random = world[random_points, ]

## although coordinates are longitude/latitude, st_intersects assumes that they are planar

nrow(world_random)

## [1] 4

#> [1] 4
random_joined = st_join(random_points, world["name_long"])

## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

plot(st_geometry(cycle_hire), col = "blue")
plot(st_geometry(cycle_hire_osm), add = TRUE, pch = 3, col = "red")

any(st_touches(cycle_hire, cycle_hire_osm, sparse = FALSE))

## although coordinates are longitude/latitude, st_touches assumes that they are planar

## [1] FALSE

#> [1] FALSE

cycle_hire_P = st_transform(cycle_hire, 27700)
cycle_hire_osm_P = st_transform(cycle_hire_osm, 27700)
sel = st_is_within_distance(cycle_hire_P, cycle_hire_osm_P, dist = 20)
summary(lengths(sel) > 0)

##    Mode   FALSE    TRUE 
## logical     304     438

z = st_join(cycle_hire_P, cycle_hire_osm_P,
            join = st_is_within_distance, dist = 20)
nrow(cycle_hire)

## [1] 742

#> [1] 742
nrow(z)

## [1] 762

#> [1] 762

plot(cycle_hire_osm["capacity"])

plot(z["capacity"])

nz_avheight = aggregate(x = nz_height, by = nz, FUN = mean)

plot(nz)

nz_heighest = nz_height %>% top_n(n = 1, wt = elevation)
canterbury_centroid = st_centroid(canterbury)

## Warning in st_centroid.sf(canterbury): st_centroid assumes attributes are
## constant over geometries of x

st_distance(nz_heighest, canterbury_centroid)

## Units: [m]
##        [,1]
## [1,] 115540

id = cellFromXY(elev, xy = c(0.1, 0.1))  #座標でraster画像のIDを取得
elev[id]

## [1] 16

raster::extract(elev, data.frame(x = 0.1, y = 0.1))  #上と同じ

## [1] 16

clip = raster(xmn = 0.9, xmx = 1.8, ymn = -0.45, ymx = 0.45,
              res = 0.3, vals = rep(1, 9))
elev[clip]

## [1] 18 24

# we can also use extract
# extract(elev, extent(clip))

# create raster mask
rmask = elev 
values(rmask) = sample(c(NA, TRUE), 36, replace = TRUE)
plot(rmask)

# spatial subsetting
plot(elev[rmask, drop = FALSE])           # with [ operator

plot(mask(elev, rmask))                   # with mask()
plot(overlay(elev, rmask, fun = "max"))   # with overlay

elev + elev

## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 72  (min, max)

elev^2

## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1296  (min, max)

log(elev)

## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : -Inf, 3.583519  (min, max)

elev > 5

## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)

library(sf)
library(raster)
library(dplyr)
library(spData)
library(spDataLarge)

## 
## Attaching package: 'spDataLarge'

## The following object is masked _by_ '.GlobalEnv':
## 
##     random_points

plot(us_states)

regions = aggregate(x = us_states[, "total_pop_15"], by = list(us_states$REGION),
                    FUN = sum, na.rm = TRUE)
regions2 = us_states %>% group_by(REGION) %>%
  summarize(pop = sum(total_pop_15, na.rm = TRUE))

## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar

us_west = us_states[us_states$REGION == "West", ]
us_west_union = st_union(us_west)

## although coordinates are longitude/latitude, st_union assumes that they are planar

plot(us_west)

plot(us_west_union)

texas = us_states[us_states$NAME == "Texas", ]
texas_union = st_union(us_west_union, texas)

## although coordinates are longitude/latitude, st_union assumes that they are planar

plot(texas)

plot(texas_union)

data("elev", package = "spData")
clip = raster(xmn = 0.9, xmx = 1.8, ymn = -0.45, ymx = 0.45,
              res = 0.3, vals = rep(1, 9))
elev[clip, drop = FALSE]

## class      : RasterLayer 
## dimensions : 2, 1, 2  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : 1, 1.5, -0.5, 0.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 18, 24  (min, max)

plot(elev[clip, drop = FALSE])

data("dem", package = "spDataLarge")
plot(dem)

dem_agg = aggregate(dem, fact = 5, fun = mean)
plot(dem_agg)

dem_disagg = disaggregate(dem_agg, fact = 5, method = "bilinear")
plot(dem_disagg)

identical(dem, dem_disagg)

## [1] FALSE

srtm = raster(system.file("raster/srtm.tif", package = "spDataLarge"))
zion = st_read(system.file("vector/zion.gpkg", package = "spDataLarge"))

## Reading layer `zion' from data source `/usr/local/lib/R/site-library/spDataLarge/vector/zion.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 11 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 302903.1 ymin: 4112244 xmax: 334735.5 ymax: 4153087
## Projected CRS: UTM Zone 12, Northern Hemisphere

zion = st_transform(zion, projection(srtm))

plot(srtm)

plot(zion)

## Warning: plotting the first 10 out of 11 attributes; use max.plot = 11 to plot
## all

#cropでくり抜き、maskで作成
srtm_cropped = crop(srtm, zion)
plot(srtm_cropped)

srtm_masked = mask(srtm, zion)
plot(srtm_masked)

zion_nlcd = raster::extract(nlcd, zion, df = TRUE, factors = TRUE) 

## Warning in .local(x, y, ...): Transforming SpatialPolygons to the CRS of the
## Raster

## Warning in spTransform(xSP, CRSobj, ...): NULL target CRS comment, falling back
## to PROJ string

plot(zion_nlcd)

dplyr::select(zion_nlcd, ID, levels) %>% 
  tidyr::gather(key, value, -ID) %>%
  group_by(ID, key, value) %>%
  tally() %>% 
  tidyr::spread(value, n, fill = 0)

## # A tibble: 1 x 9
## # Groups:   ID, key [1]
##      ID key    Barren Cultivated Developed Forest Herbaceous Shrubland Wetlands
##   <dbl> <chr>   <dbl>      <dbl>     <dbl>  <dbl>      <dbl>     <dbl>    <dbl>
## 1     1 levels  98285         62      4205 298299        235    203701      679

cycle_hire_osm_projected = st_transform(cycle_hire_osm, 27700)
raster_template = raster(extent(cycle_hire_osm_projected), resolution = 1000,
                         crs = st_crs(cycle_hire_osm_projected)$proj4string)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in Proj4
## definition

plot(cycle_hire_osm_projected)

ch_raster1 = rasterize(cycle_hire_osm_projected, raster_template, field = 1)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in Proj4
## definition

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum OSGB 1936 in Proj4 definition

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)

##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.